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We study the evolution of diffuse elastodynamic spectral energy density under the influence of 
weak nonlinearity. It is shown that the rate of change of this quantity is given by a convolution 
of the linear energy at two frequencies. Quantitative estimates are given for sample aluminum and 
fused silica blocks of experimental interest. 



I. INTRODUCTION 

Weak nonlinearity is known to cause a redistribution of 
elastodynamic spectral energy density. Energy present at 
one or two base frequencies will migrate, under the influ- 
ence of nonlinearity, to higher harmonics and to sum and 
difference frequencies. If those sum and difference fre- 
quencies originally have little or no energy, their presence 
can be a sensitive measure of nonlinearity. The method 
has been shown to be capable of detecting flaws in oth- 
erwise linear-response specimens. It has been applied to 
crack detection [l| and to assessing the quality of adhe- 
sive bonds jjj. Much of this work has been qualitative, 
with emphases on how the nonlinear effects scale with 
field amplitude. The geometry is such, and the fields 
sufficiently uncontrolled, that attempts to quantify the 
nonlinearity in absolute terms have not been feasible. 

In special circumstances, with high amplitude and long 
distance plane wave propagation, it is possible to moni- 
tor harmonic generation, and thereby assess nonlinearity, 
and do so quantitatively. It has been suggested that non- 
linear Rayleigh wave propagation may be a useful NDE 
technique to measure surface properties Q . Lamb waves 
were also suggested for this purpose 0|. 

As an alternative to plane wave configurations diffuse 
elastodynamic fields might be appropriate for NDE mea- 
surements under insufficiently controlled conditions. Dif- 
fuse fields span a broad range of applications from op- 
tics and microwaves to ultrasonics (see Q and references 
therein). Gross properties of a specimen can be evalu- 
ated from integral field parameters @. Extension of dif- 
fuse field theory to include nonlinear effects is receiving 
increasing experimental Q and theoretical Q attention. 
Its application to evaluation of material properties may 
be a future research focus. 

The purpose of this paper is to describe effects of weak 
nonlinearity on the energy spectrum of a diffuse field, as 
one of the integral quantities available for experimental 
measurement. With a view toward eventual application 
in the context of NDE, but with an interest in the general 
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problem, we take elastodynamics in a nonlinear solid as 
an underlying physical system for our diffuse field model. 

We start with the governing equations for the chosen 
field type in Section II, and develop them into a system of 
ODEs with quadratic and cubic nonlinear terms describ- 
ing time-evolution of the modal amplitudes. In Section 
III we discuss statistical assumptions for the normal fre- 
quencies and modes of the solid. In Section IV the non- 
linear equations are treated by means of regular pertur- 
bation theory and averaged to obtain the spectral power 
transfered into a frequency band. Results independent 
of the underlying physical nature of the nonlinearity are 
considered. In Section V we specialize to the case of an 
isotropic homogeneous elastic body and an initial field 
generated by two transient narrow-band signals centered 
at different frequencies. Quantitative estimates for the 
nonlinear energy transfer into double and combination 
frequencies are provided and discussed. Conclusions are 
presented in Section VI. 



II. GOVERNING EQUATIONS 

Suppose we have a finite hyperelastic body occupying 
volume Vq in its natural state referenced by coordinate 
x, and having material density po (x). Strain energy den- 
sity of the body for a displacement field u with a corre- 
sponding Green's tensor Eij — 1/2 (mj + Ujj + Uk,iUk,j) 
is described to leading orders in strain by 

W=^C ijkl (yL)EijE kl 

■f" ~^ Dij]~i ran (x) E k i E 7nn 

~t~ ~^Eij k l rnn pq (x) E{j E k lE mn Epq -f- . . . , (1) 

where C, D and F are linear, second and third-order 
nonlinear elastic tensors respectively. 

Let the body be subjected to external forces and trac- 
tions that stop acting after a certain cutoff time to = af- 
ter which the field in the body freely evolves without dis- 
sipation under zero-displacement boundary conditions. 
According to the Hamilton's Principle, the weak form 
of the governing equation for the field evolution can be 
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written as 



llvo 



(—(jjiSui j — poiiiSui) dx 3 dt = 0, (2) 



where a is the first Piola-Kirchoff stress tensor. It is 
related to the strain energy density Q via deformation 
gradient Fij — Sij + mj 



a 3i (x,u) = dW/dF h 



(3) 



To distinguish parts of equation © responsible for lin- 
ear and leading-order nonlinear behavior of the system 
we expand the stress tensor a in powers of ||u|| up to the 
third-order, and label the respective linear and nonlinear 
operators as L and N. Equation is then written as 



Sui [poUt - Liu) + SuijNijU 



dx 3 dt = 0. (4) 



A complete set of eigenvalues cj„ and their normalized 
eigenfunctions u" is considered to be known for the linear 
operator as solution of 



L,-u" 



_d_ 



9 n 



(5) 



For clarity of notation we introduce composite Greek 
indices with implied summation upon them: 



gij (x) hij (x)dx 3 , 



v 



and denote the first partial derivatives of the field as sep- 
arate functions: it Q ={ x .i.j} = dui (x) /dxj. 

The displacement field allows decomposition upon the 
eigenfunctions with time-dependent modal amplitudes dk 



cikc] 
.10]: 



(x, <)=$>* (*)«*(*)• 



By employing this representation in equation and 
using the eigenfunctions as a set of test functions Su, we 
restate the governing equation as 



d k + uj k d k 



/"] N a U^p^ dk di d„ 



^ N a/ 3 7 sU^gdkdid m d n , (6) 



J,' 



where Uj^fi 



jkhnn 



u~upu-, and U^ s 



i k a u l pu™u%. All 

the specifics of particular type of nonlinear behavior are 
now contained in modal matrices NU. 

Only symmetrical with respect to the last two indices 
parts of matrix A r Q , ( g 7 and symmetrical with respect to the 
last three indices parts of matrix N a p 7 s survive modal 
summation in 0. Without loss of generality we put 



these matrices equal to their corresponding symmetrical 
parts, and write them in Greek index notation as 



N. 



i 
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NijklmnS (x - x') S (x - x") . 



and 



A' ck— {x,2 J< 7'}/3={x / ,fc,Z}7={x // ,m,n}(5={x /// ,p,q} 

= ^NijklmnpgS (x - x') S (x - x") 5 (x - x'") . 

Nonlinear terms of the stress (j^J yield directional tensors 
in the above formulae: 

A' ijklmn 

— Cijln^km C j n kl&im C * jlmn^ik ~^ Dijklmm (7) 

and 

A ijklmnpq 

— Cjlnq3ik&mp ~t~ Cj n lqSi m 8kp ~\~ Cjq n lSip5 7n k 
~t~ D jlmnpq&ik ~t~ ^jnklpq^im ~t~ ^jqmnkl^ip 
~t~ DijlnpqSmk ~t~ Dijlq mn Spk ~\~ Dijqnkl&pm 
~\~ Fijklmnpq- 

Due to the major and minor symmetries of elastic tensors 
C, D and F [9] both matrices N can be identified as fully 
symmetric: 



(8) 



A quantity which we term linear energy stored in a 
single mode k at time t is 



Ek= 1 -(d 



1 1 "14 



(9) 



In the absence of nonlinearity it is equal to the total 
energy of a mode, and is constant over time. It can be 
written in terms of the complex amplitudes ipk- 



E. 



(0) 



1 



(10) 



The amplitudes ipk arise from the action of external forces 
and tractions prior to the cutoff time t . They describe 
evolution of the linear part of the field by means of the 
modal amplitudes d k °^ that are found as solutions of the 
linearized version of equation ©, with matrices N put 
to zero [ill ]: 



4 0) = a k (t) = ^ k e~ Wkt . 



(11) 



The energy flow due to nonlinearity, E k — Iljt , is 
obtained from the governing equation 10, with modal 
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power input being 

life = — iV Q/ 3 7 ^2 Uaf3^dkdld m 
in .1 

- N a p lS U^d k did m d n . (12) 

m,l,n 

It must be noted that the energy quantity E k © is not 
conserved, because it is not a true energy in the non- 
linear case. Analysis of the strain energy density re- 
veals that additional terms — A r Q , ( a 7 ^2 m t Uj^dkdidm/3— 

N aP78Y, m ,i, n U af^8 d k did m d n / '4 must be added to E k to 
produce a quantity that is conserved. However, the terms 
involve summation upon modal amplitudes other than 
that of the mode at hand, and do not allow simple inter- 
pretation in terms of a single mode. Thus they are not 
used. 



III. STATISTICAL MODEL 

We restrict ourselves to a class of systems for which 
the field excited by external forces and tractions has a 
fully diffuse nature. In experiments such fields are prac- 
tically realized, for example, in an elastic solid of a classi- 
cally chaotic shape, and have statistical properties close 
to or indistinguishable from those of a field described by 
a random Hamiltonian (l2llT^ |. The normal modes of the 
system are taken to be centered Gaussian vectors with a 
certain spatial correlation, as was first theoretically con- 
jectured 13, and later numerically and experimentally 
verified [13 . We assume that the mean density of states 
of the normal frequencies is given by the function D (u>) 
in the form of Weyl-series [12], and frequency-frequency 
correlations can be neglected. 

Since experimental identification of a particular mode 
is complicated by the "missing level" effect or modal over- 
lap we choose to pursue calculation of average spectral 
density rather than individual modal amplitudes. In or- 
der not to distinguish between individual modes at the 
frequencies of interest, and thus deal with the average 
quantities, we limit observation time of the system i to 
be less than the corresponding break (Heisenberg) time, 
iff = 2-kD. On the other hand, time t is considered larger 
than transition (ballistic) time in the solid (mean time 
between two successive scattering events at the bound- 
aries), ti = l/c, so that diffuse regime of the field is estab- 
lished. I and c stand for characteristic diameter of and 
wavespeed in the solid. Putting together the two bounds 
yields t; <C t <C iff, a condition that can be experimen- 
tally realized. 

In the framework of the adopted statistical model the 
modes u k of the linear operator © are centered Gaussian 
random vectors with variance given by 

(u* (x) m™ (x')> = SknKij (wfe, x, x') , 

where (...) represents ensemble average, and K is a 
smooth frequency-dependent correlation matrix. Pair- 



wise correlation of the first partial derivatives of the 
modes 

(uluf) = S kn K' aP K) (13) 

is readily obtained from 

^C={x,M}/3={x'j\m} = 9 2 Kij (cj,x,x') / dxidx' m . (14) 

Orthonormality of the modes imposes a normalization 
condition upon the correlation matrix: 

/ p (x)K H (u,x,x)dx 3 = 1. (15) 

JVo 

From the mode statistics the complex amplitudes 4>k 
describing initial field are found to be centered Gaussian 
random numbers with variance 

(ipklpn) = 0, -D (u! k ) (ipllpn) = S nk £ (Wfe) . 

The amplitudes ip k relate to the linear energy iflOl) . and 

give the smooth function e (uj k ) — D (oj k ) (e^ / wjA an 

interpretation as a frequency-normalized spectral energy 
density. The linear amplitudes a k 111 111 form a centered 
Gaussian ensemble as well. Their pairwise time correla- 
tion is calculated as 

D (u k ) (a k (i) a n (r)) = 8 kn e (uj k ) coscjfc (t - r) . (16) 

Both, mean density of states and correlation ma- 
trix, can be calculated in terms of the average Green's 
function: D (u) = 2u)% (Tr (G /tt, and K (u>) = 

For time scales under consid- 
eration, in particular t\ <§; iff, the characteristic wave- 
length on the frequencies of interest is much smaller than 
the diameter of the solid: X/l <C 1. The leading order 
nonlinear contribution thus comes from the bulk rather 
than near-boundary region of the solid. This allows us to 
neglect effects of the latter, and approximate the exact 
Green's function G in the solid by the Green's function 
in an unbounded medium G°°. This approximation im- 
plies that K has an infinite correlation radius inherited 
from G°°, and leads to formal integral divergence in the 
calculations of the following sections. To mend the prob- 
lem we consider the scattering of the field inside the solid 
as a diffuse process with a free mean path on the order 
of the diameter I, which now provides a finite correlation 
radius for the model: (G) = (G°°) e -| x ~ x 'l /! 01- The 
ansatz is justified as the final results turn out not to de- 
pend on the specific choice of the value of I, as long as it 
stays much greater than the wavelength. 

IV. ENERGY SPECTRUM EVOLUTION 

We assume nonlinear effects to be small and seek solu- 
tion of equation by utilizing small perturbation the- 
ory, and expanding amplitudes d k in orders of magnitude 
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dk = o>k + bk + (•••) , where a k are given by (11 II . Next or- 
der amplitudes 6^ arise from the presence of nonlinearity, 
and are determined solely by the linear field: 

b k a) = / V in " fe( *- r) 

Jtn U k 



E^ 7 ^™a ; (r)a m (r) 

■ m,i 

■ 53 N «0nS u *(hS a i ( T ) a »> ( r ) a « (t) 

m,l,n 



To obtain successive corrections to the average spectral 
density of the power flow we expand l(T2|l in a series of 
magnitudes of a k : 

n («*, t) = z? (w fc ) (m (t)) = 53 n< n ' (a;*, t) . (17) 



As mentioned in Section II, since the power flow arises 
from the nonlinear mode coupling, the modal energy of 
the linear field Ijlflll is conserved: II^ ) (u,t) — 0. For the 
assumed statistics of the amplitudes a k (IKifl the first- 
order correction to the average power density is zero as 
well: 



n (1) ( Wfcj f) 

=D (uj k ) J] N aPl (U$%) (d fc (t) ai (t) a m (t)) = 

7 



= 0. 



1,171 



The power flow expansion itTTjl starts with II' 2 ^ as the 
leading term. We express the fourth moments of the 
amplitudes a k as double products of their pairwise cor- 
relations iflfil) , and obtain the power flow in terms of the 
energy density: 



n( 2 ) (w fc) t) = D K) 53 N a01 N um (U^U k ™ v m ) e' ( Wm ) e' (u m ) 



smuj k t 



J 53 N a01 N vm (UffiU?™) e> (un) 



e> (co m ) ± 2e' ( Wfc ) ^ 



sin (w fe - w ; ± w m ) f sin(w fe + cj z ±w m )* 



OJfc — LJ; ± ijJ T , 



UJk +UJi±UJr, 



(18) 



where e' (u) — e(ui)/D(u), is the average normalized 
energy of a single mode. The modal sum is evaluated as 
a frequency integral with integrand weighed by the mean 
density of states: 

r+oo 

53 fn = / / D i^n) duj n 

We note that matrix N a pys responsible for the cubic non- 
linearity does not enter l|T8jl . for its contribution is pro- 
portional to (dfc (t) ai (t) a m (t) a n (t)) = 0. 

The centered Gaussian statistics of the modes allows 
averages in the form of (u^^U^^J to be expressed as 

triple products of the pairwise correlations 111 -it . Com- 
bining terms with the same modal indices we formally 
write (with no implied summation on Latin indices) 

N a (3~fN vm ( U a ^ U v ™) 

+ N 2 (uji,uj m )8 lk +N 3 (uj k )6 k iS km . (19) 

The coupling functions N are defined later. 

Expression Ill8|l with the factors l|19fl inserted, though 
bulky and cumbersome to analyze, gives the leading term 
in the power flow due to a weak nonlinearity. In exper- 
imental practice, however, it is not uncommon to deal 



with the fields that have the spectral energy density e 
varying smoothly on frequency scales Aw greater than 
the averaging bandwidth of the limited observation time: 
Aujt 3> 1. In this case ijlSjl simplifies further. All the 
terms directly proportional to rapidly oscillating sine fac- 
tors yield negligible averages: 

sin (uj n t) /u} n = D (w„) O (1/Awt) , 

and the modal sums involving such factors are evaluated 
as follows: 

53 fn sin [K -u)t]/ (u) n - uj) 

n 

=ttD (w) / (u) [6 (w) + O (l/Aut)] , 

where 9 is a unit step function. Another simplification 
arises if we pay attention only to the frequencies that car- 
ried no energy initially: e (u) — 0. For such frequencies 
a weak change in the energy density due to II^ 2 ^ is not 
masked by a strong initial linear field, and is convenient 
for experimental measurement. 

With the simplifications mentioned the power flow ex- 
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pression reduces to 

_ r+oo 

n( 2 )( w ,t) = -^ / D(u>)X (lj,u',\w±w'\) 
z ± Jo 

x e(\oj± w'\)e(w')du'. (20) 

The only remaining coupling function is given by con- 
traction of the nonlinear and correlation matrices: 

= N a0J N vm K' av A*^ (w 2 ) if' 7 „ (073) . (21) 
We note that 

N (u,u', \oj + u>'\) = N (oj, w", |w - w"|)| w „ =w+w ,. 

The power flow in the form of <2(H allows simple inter- 
pretation: the energy transfered into a given frequency w 
comes from all pairs of frequencies 10' and to" that have 
the given frequency as a combination, i.e. equal to their 
sum or difference: uj — \u>'±u>"\. The qualitative re- 
sult is in agreement with and could have been expected 
from an elementary theory of nonlinear oscillations [Tfj |. 
The symmetry of the nonlinear matrix (JBJ leads to the 
symmetry of No with respect to any interchange of its ar- 
guments. The coupling strength of any triad of frequen- 
cies is thus independent of the energy transfer direction. 
Nevertheless, the power itself exhibits an overall trend of 
the energy to be transfered up the frequency spectrum, 
as it is proportional to the density of states at the target 
frequency D (u>). 

The expression 112011 implies energy growth that is, if 
e's are approximately constant in time, proportional to 
elapsed time t. In discrete spectrum systems such be- 
havior is found when a triad of frequencies is locked in 
internal resonance, producing secular terms in the solu- 
tion obtained by regular perturbation theory. This is the 
behavior observed in our case because of the finite time 
i « (g: since individual modes are not resolved, any 
combination frequency produced by the source frequen- 
cies in the frequency range of interest is indistinguishable 
from at least one of the normal frequencies of the system, 
which hence lies in effective internal resonance with them. 



Appendix A): 



^■(Ax, w ) = ^ e -l Ax l/ ! 



1 

AT 



1/cf + 2/c? 



3 



ijo (*J |Ax|) + 4io (*t|Ax|) 
c i c t 

{Sij/3 — AxiAxj) 



iia (fc/|Ax|)-ij 2 (fc t |Ax|) 



(22) 



where Ax = x — x' is a separation variable, Ax = 
Ax/|Ax|, M is total mass of the solid, cj and c* are 
the longitudinal and transverse wavespeeds respectively, 
and j n is the spherical Bessel function of order n. The 
mean density of states for clamped boundary conditions 
is calculated in [lflj : 

D(co) =^ [l/cf + 2/c 3 ] 



So , 
8ncf 



2 + {ci/c t y +3( Cl /c t y 
(d/ctf + l] +0(l/c), 



with Sq being the surface area of the solid. 

For the model of nonlinearity at hand appropriate for 
an isotropic solid, described by the five-constant theory 
[2C| . and correlation matrix K provided by ll22|l . the cou- 
pling function No t2H assumes the form 



No {lo,lu',\uj~J\) 

7T C 

~ ( 7; 3 + 2 7t 3 ) 3 V M ' 



1 

UJ 



— N (uj'/oj) . 



N is a dimensionless function of the frequency ratio and 
material properties that characterize nonlinear coupling 
strength between the source and target frequencies. The 
function is independent of the linear dimensions of the 
body and possesses the symmetries according to (JBJ: 
N(u'/u) = N(w/u') = N(l- uj'/uj). An analytical ex- 
pression for N is available, but too bulky to be presented 
here. Details of the calculations are found in Appendix 
A. 

As a sample distribution of the (linear) energy density 
we take two Gaussian peaks with half- width Alo and total 
energies centered at frequencies u>\ and u>2'- 



V. ESTIMATES FOR ISOTROPIC 
HOMOGENEOUS ELASTIC SOLID 



»={1,2} 



V2ttAw 2 



= -(w-tj i ) 2 /2Ati 



(23) 



We return to elastodynamic displacement fields, and 
specialize to the case of an elastic body composed of 
known isotropic homogeneous material, the case that 
holds premium experimental and theoretical interest. 

The linear part of the Green's function in the un- 
bounded medium yields the correlation matrix K (see 



The width of the peaks is restricted by condition Aut ^> 
1, as imposed by applicability requirements for the sim- 
plified power input expression i2()\) . It is additionally 
assumed that the peaks do not overlap: \<u)\ — w 2 | 3> Aw. 
The suggested form of the spectrum closely models the 
combined spectra of two narrow-band signals that might 
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Figure 1: Dimensionless coupling function in aluminum 
(thick line) and fused silica (thin line). 

be used in an experiment. By adjusting the carrying fre- 
quencies the features of nonlinear mode coupling can be 
investigated, and the form of N measured. 

As a parameter suitable for characterization of the non- 
linear energy transfer strength we consider transfer times 
defined as the formal time required for entire mean en- 
ergy of the source peaks 12311 to be transfered into com- 
bination frequencies, provided small perturbation the- 
ory and equation (1201 holds: t± t 2 — -Ei,2 / and 

ti±2 = \fE\Ei I . Here n( 2 ) represents the total 

power input into frequency band supporting the result- 
ing peak at the combination frequency. According to the 
definition ti t 2, ii+2, and ii_ 2 characterize energy transfer 
into double, sum and difference frequency respectively. 
Calculations yield these times to be 

. = l gigg J_ 

7r Mc 2 D Quji ± w 2 |)No (\u>i ±u> 2 \ ,a;i,a;2) mj ±2 ' 

(24) 

t 2= 2 <4a i_ 

7rMc 2 D(2cji :2 )N (2u;i, 2 , ^1,2, ^1,2) mf 2 

As a measure of absolute strength of the linear field 
we choose squares of the Mach numbers: m\± 2 = 
^/E\E^/Mc 2 and m\ 2 = E lt2 /Mc 2 . As seen from JMJ 
stronger nonlinear effects (shorter times) require higher 
Mach numbers, i.e. stronger linear fields. 

For numerical estimates we take an aluminum block 
and a fused silica block with the mechanical properties 
listed in Tabled and choose the transverse wavespeed for 
each material as the characteristic one: c = q. Due to 
the inversion symmetry we evaluate the coupling function 
N for frequency range < uj'/uj < 1 only, and plot it in 
Figure Q Immediately a significant feature of the plot 
comes into view: the sharp peaks and discontinuities of 
the coupling function at frequency ratios of 0.25 and 0.75 
for aluminum, and 0.17 and 0.83 for fused silica. 

The feature may be understood by considering nonlin- 
ear combinations of plane waves. We recall that from a 



plane wave perspective appropriate in the limit X/l <C 1, 
in addition to internal frequency resonance between the 
source frequencies u>' and u>" = \u> — u>'\ and the tar- 
get frequency u, there is an additional requirement for 
wavevector resonance. We note that there are two plane 
wave types in the system, namely longitudinal and trans- 
verse. If the types of the three waves (at u), u/ and oj") are 
identical, wavevector resonance (k(ui") — k(w) + k(u/) 
with |k| = lo/c) is always possible, and demands that the 
three wavevectors are parallel or antiparallel. If the wave 
type for one of the frequencies is different from the other 
two, then wavevector resonance is not always possible; it 
depends on the frequency ratio. If it is possible, there 
will be a nontrivial angle between the wavevectors. The 
transition between possibility and impossibility occurs at 
certain special frequency ratios, lo/lo' — (l±c t /q)/2. 
At these ratios wavenumber resonance occurs with par- 
allel or antiparallel wavevectors. The special frequency 
ratio depends solely on the Poisson ratio of the mate- 
rial. For non-dispersive single wavespeed systems, when 
q = c t , the ratios become and 1. In this case the effect 
is not observed, as one of the source frequencies needs 
to be zero [24|]- However, for certain dispersive single 
wavespeed systems the peaks and discontinuities might 
still be found, if the dispersion equation of the system 
is such that the wavenumber resonance is possible. It is 
also worth noting that the discontinuity position at lead- 
ing order is independent of the field strength. Thus field 
calibration is not required in experimental measurements 
in order to observe this characteristic feature, for only 
relative values of N are needed. This makes the method 
convenient to use in the light of its eventual application 
to NDE. 

In order to see whether the presented theory is appli- 
cable for experimental verification and eventual use, we 
need to provide numerical estimates for transfer times 
(I24T) . For this purpose we choose the solid to be equiva- 
lent in volume and surface area to a cube with a side of 
7cm. The typical values of the transition (ballistic) time 
in the system are of the order 10 [is. The carrying fre- 
quencies of the two narrow-band signals are taken to be 
u>i /2ir = 500 kHz and w 2 /27r = 600 kHz, so that the break 
(Heisenberg) time is of the same order as the absorption 
times of 100 ms common to experiment. The transfer 
times are calculated for the aluminum block and plot- 
ted in Figure El as functions of the linear field strength, 
characterized by Mach number. 

It is feasible to generate linear diffuse fields with RMS 
elastic strains corresponding to Mach numbers of order 
10 -5 that yield transfer times of order 10 3 s (see Figure 
|2j. Thus the observation times of the order of the ab- 
sorption time yield energy densities at combination fre- 
quencies some 10 4 times weaker than the source signal 
energy, i.e. RMS strains of 10~ 7 . These are easily de- 
tectable. We also note that the strain ratio (combination 
field to initial field) is much smaller than unity, so the 
use of regular perturbation theory is valid for the given 
time scale. 
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Table I: Mechanical properties. 





Po (kg/m 3 ) 


c; (m/s) 


Ct (m/s) 


A (GPa) 


At (GPa) 


A (GPa) 


B (GPa) 


C (GPa) 


Aluminum" 


2720 


6100 


3090 


49.1 


26.0 


-320 


-198 


-190 


Fused Silica 


2200 6 


5700 


3750 


9.60 6 


30.9 6 


-44 c 


93 c 


27 c 



"Alloy D54S, Smith et al. [2l| 
6 Drumheller E3l 
c Bechmann et al. l23ll 



t(s) 




1 1 m 1±2 ,mi 2 

lxl0~ 5 3xl0~ 5 7xl0~ 5 



Figure 2: Transfer times in aluminum block for fixed o;i/27r = 
500 kHz and uj 2 /2tv = 600 kHz. 



t(s) 




Figure 3: Transfer times in aluminum (thick lines) and fused 
silica (thin lines) block for varying source frequency ratio, and 
fixed U2/2ir — 550 kHz and mi±2 = 10 . 

Finally, we calculate transfer times in aluminum and 
fused silica as a function of source frequency ratio for a 
fixed Mach number (see Figure EJ ■ The frequency de- 
pendence of the nonlinear coupling strength in terms of 
the transfer times reveals the frequency dependence of N. 
The dips in the transfer times correspond to the peaks 
of the coupling function, and might be used to estimate 
Poisson's ratio of the material, or as a signature of nonlin- 
earity. We notice that the energy transfer into near-zero 
frequencies is inefficient, as manifested by high difference 
frequency transfer times in the vicinity of u>i/lo2 = 1. 



VI. CONCLUSION 

In the present work we have obtained a formula for the 
evolution of the average linear energy spectrum (l20l of 
a weakly nonlinear system. Except for definition of the 
modal linear energy Ek (©, and statistical properties of 
the normal frequencies luj. and modes u fc (JjJ, the result 
has no reference to the physical nature of the dynamic 
system governed by nonlinear equations © • With proper 
generalization of the matrices NU describing the physical 
nature of the nonlinearity it might be applicable beyond 
the scope of elasticity. 

We observe that the energy redistribution occurs for 
triads of frequencies, with one being the sum or difference 
of the other two. The average power input into a narrow 
frequency band is found to be a cumulative effect coming 
from such interactions in the initial energy spectrum of 
the field, and to be proportional to the convolution of the 
energies stored at two frequencies. Relative weight of the 
interactions is given by the frequency-dependent coupling 
function No (12 H , and is calculated as a contraction of the 
nonlinear and correlation matrices N and K'. 

The case of chief experimental interest involving an 
isotropic homogeneous elastic body with nonlinearity 
given by the five-constant theory and an initial energy 
spectrum consisting of two narrow-band signals is dis- 
cussed in detail. We find that characteristic times for 
the full energy transfer from the source into combina- 
tion frequencies depend on the ratio of the source fre- 
quencies and exhibit characteristic dips and peaks and 
discontinuities at special frequency ratios corresponding 
to wavenumber resonance. The position of the resonance 
depends on Poisson's ratio. We also find that the transfer 
times are reciprocal to the strength of the initial signal 
given by the square of the corresponding Mach number. 

The current theory is derived for zero-displacement 
boundary conditions of a clamped solid, and is imprac- 
tical for experimental realization. However, for the fre- 
quencies of interest the main contribution to the mix- 
ing comes from the bulk of the solid. Since the near- 
boundary regions play a lesser role, the authors believe 
that the energy spectrum evolution for the traction-free 
elastic solid would exhibit behavior similar to that dis- 
cussed here, and thus be accessible for experimental ver- 
ification. 



8 



[1] V.Yu. Zaitsev, A.M. Sutin, I.Yu. Belyaeva, and V.E. 
Nazarov, J. Vibration and Control 1(3), 335 (1995); 
Alexander M. Sutin and Dimitri M. Donskoy, in Proceed- 
ings of the 8th International Symposium on Nondestruc- 
tive Characterization of Materials, Boulder, Colorado, 
1997, edited by Robert E. Green Jr. (Plenum Press, New 
York, 1998) 

[2] A. Wegner, A. Koka, K. Janser, U. Netzelmann, S. 
Hirsekorn, and W. Arnold, Ultrasonics 38, 316 (2000); 
M. Rothenfusser, M. Mayr, and J. Baumann, ibid. 38, 
322 (2000) 

[3] D.C. Hurley, J. Acoust. Soc. Am. 106(4), 1782 (1999) 
[4] W.J.N, de Lima and M.F. Hamilton, J. Sound Vib. 

265(4), 819 (2003) 
[5] Hans-Jiirgen Stockmann, Quantum Chaos: an Introduc- 
tion, (Cambridge University Press, 1999), Chap. 2, 6 
[6] D.J. Pine, D.A. Weitz, P.M. Chaikin, and E. Her- 

bolzheimer, Phys. Rev. Lett. 60(12), 1134 (1988); 

Richard L. Weaver and Wolfgang Sachse, J. Acoust. Soc. 

Am. 97(4), 2094 (1995) 
[7] Johannes F. deBoer, Ad Lagendijk, Rudolf Sprik, and 

Shechao Feng, Phys. Rev. Lett. 71(24), 3947 (1993) 
[8] R. Bressoux and R. Maynard, Europhysics Letters 50(4), 

460 (2000); S.E. Skipetrov and R. Maynard, Phys. Rev. 

Lett. 85(4), 736 (2000); S.E. Skipetrov, Phys. Rev. E 63, 

056614 (2001) 

[9] RW. Ogden, Nonlinear Elastic Deformations (Ellis Hor- 
wood Ltd., Chichester, 1984), Chap. 3.4, 4.3, 6.1 
[10] Under zero-displacement boundary conditions solution of 
the nonlinear problem belongs to the function space with 
a basis formed by the eigenfunctions of the linear opera- 
tor, thus allowing the proposed decomposition. For zero- 
traction boundary conditions such decomposition is no 
longer valid, since the eigenfunctions no longer span the 
entire solution space of the nonlinear problem. 
[11] R.L. Weaver, J. Sound Vib. 94(3), 319 (1984) 
[12] R.L. Weaver, J. Acoust. Soc. Am. 85(3), 1005 (1989) 
[13] C. Ellegaard, T. Guhr, K. Lindemann, H.Q. Lorensen, 
J. Nygard, and M. Oxborrow, Phys. Rev. Lett. 75(8), 
1546 (1995); C. Ellegaard, T. Guhr, K. Lindemann, J. 
Nygard, and M. Oxborrow, ibid. 77(24), 4918 (1996); A. 
Andersen, C. Ellegaard, A. D. Jackson, and K. Schaadt, 
Phys. Rev. E 63, 066204 (2001) 
[14] M.V. Berry, J. Phys A 10(12), 2083 (1977); P. O'Connor, 
J. Gehlen, and E.J. Heller, Phys. Rev. Lett. 58(13), 1296 
(1987) 

[15] Steven W. McDonald and Allan N. Kaufmann, Phys. 
Rev. A 37(8), 3067 (1988); K. Schaadt, T. Guhr, C. 
Ellegaard, and M. Oxborrow, Phys. Rev. E 68, 036205 
(2003) 

[16] The expressions are obtained by calculat- 
ing ensemble average of the imaginary part 
of the Green's function 3Gy (x,x',w) = 

* £»=i "i* ( x ) «j° ( x ') 5 (<" - w «) E.N Economou, 

Green's functions in quantum physics (Springer- Verlag, 

New York, 1979), Chap. 1 
[17] V.N. Prigodin, B.L. Altshuler, K.B. Efetov, and S. Iida, 

Phys. Rev. Lett. 72(4), 546 (1994); V.N. Prigodin, ibid. 

74(9), 1566 (1995) 
[18] Leonard Mierovitch, Elements of Vibration Analysis 

(McGraw-Hill, New York, 1986), Chap. 10 



[19] Marc Dupuis, Robert Mazo, and Lars Onsager, J. Chem. 
Phys. 33(5), 1452 (1960) 

[20] F.D. Murnaghan, Finite Deformations of an Elastic Solid 
(Chapman & Hall, New York, 1951), Chap. 4 

[21] R.T Smith, R. Stern, and RW.B Stephens, J. Acoust. 
Soc. Am. 40(5), 1002 (1966) 

[22] D.S. Drumheller, Introduction to Wave Propagation in 
Nonlinear Fluids and Solids (University Press, Cam- 
bridge, 1998), App. B 

[23] R. Bechmann and R.F.S. Hearmon, in Elastic, Piezo- 
electric, Piezooptic, Electrooptic Constants, and Nonlin- 
ear Dielectric Susceptibilities of Crystals, edited by K.H. 
Hellwege and A.M. Hellwege, Landolt-Bornstein, New 
Series (Springer- Verlag, New York, 1969), Vol. HI/2 

[24] For this reason a discussion of weak nonlinearity in such 
simple systems would offer little insight and few rich be- 
haviours. 

[25] R.L. Weaver, J. Mech. Phys. Solids 38(1), 55 (1990); J.D. 

Achenbach, Wave Propagation in Elastic Solids (Elsevier, 

Amsterdam, 1999), Chap. 3 
[26] B. Eckhardt, U. D6rr, U. Kuhl, and H.-J. Stockmann, 

Europhys. Lett. 46(2), 134 (1999) 



Appendix A: CORRELATION MATRIX AND 
COUPLING FUNCTION FOR ISOTROPIC 
HOMOGENEOUS MEDIUM 

Elastic deformation of an isotropic homogeneous solid 
is described by the five-constant theory in terms of the 
Lame constants A and [x, and nonlinear coefficients A, 
B and C [13|- According to the theory the linear and 
second-order nonlinear elastic tensors have the form of 

Ciju =X5ij5 k i + 2(j,$ ijkl , 

(Al) 

Dijklmn = 2-A<&ijfci mn + 2CSijSkl5 mn 

+ 2B (<%$;Lm + ran ^mn^ijkl) ' 

The elementary isotropic tensors $ are given as follows: 



ijklmn 

= — (SikSjmSln + SikSj n Si m + SuSj m Skn + $ilO~jnO~km 
O 

+ SimSjkSin + SimSjlSkn + 0~inSjkO~lm + 8in.5jl§km) ■ 

By substituting tensors ijAlll into the strain energy defi- 
nition JQ| , particular form of directional tensors Nijkimn 
{3> is derived. 

The Green's function in the unbounded medium G°° 
is calculated from its spatial Fourier-transform [2|j: 

r«x f n . A _ PiPj , ~ PiPj A — ty/ItyI 

(Pj^J — 2 + ~^T1. 2' P ~ p/ p • 
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Direct integration of the above expression and its subse- 
quent normalization flU yields particular expression for 
correlation matrix K {22} • As two of its particular limit 
cases we note, first, the known autocorrelation function 
for scalar Helmholtz equation: jo(k|Ax|) H obtained 
by letting c = cj = Cj. And second, the autocorrelation 
function for purely transverse field - such as, for example, 
electromagnetic field j2^| - obtained by letting cj — > oo. 

The correlation matrix of the first partial derivatives 
of the modes is derived from K by means of itlljl . It is 
expressed in terms of known directional tensors H, char- 
acteristic wavenumber k = lo/c, and wavespeed ratios 
li.t = c/ci jt \ 



K' 



a— {x,z,/}/3={x / ,j,m} 



- _L p -|Ax|/i 

M 

^T9^ k2 E 7 p 5 4^a(7 P fc|Ax|). 



It + 27? 



i={0,2,4} 
P = {l,t} 



(A2) 



The sum over longitudinal and transverse wave types is 
denoted as p = {l,t}. Directional tensors H are defined 
by the following expressions: 



tt(1,0) _ - r r 
11 ijmn ~ ^OijOmn 



1 



H 



(t,0) 



11 ijmn ~ u ij^mn 



H. 



(*,2) 



tjmn 



' ijmn 1 



H 



(M) 

ijmn 



-H 



(M) 

ijmn 



-O 2 - 



o 4 



The quadruple, composite quadruple and 2 4 -order direc- 
tional moments Q moments yield zero values when inte- 
grated upon all spatial directions or contracted upon any 
two pairs of indices: 



Q% 

* I'jmn 



) 4 

J ijmn 



- Ax, A 



.1 1 



3 

1 

15 



) 2 



The coupling function 12111 is given by contraction of 
the directional tensors of the nonlinear matrix with 
correlation matrices of the first derivatives l|A2l) . It is 
used to find the dimensionless function 



N= E (7 P 7,70 4 HM c r (-7-), 

a,fa,c={0,2,4} 
p,q,x = {l,t} 

where constants H are 



TOlP^ _ AT , , AT fJ-(P.°) rr{l,b) rr(t,<=) 



and the integral / is given by 



,7p7q 



ri 11 — 01 



+00 



z 2 e -z/Z-f x k 



xj a [^nz)j b (^\i-n\ z )j c (z)dz. 



(A3) 



For a finite correlation radius of the coupling matrix I, 
the integral additionally depends on the target frequency 
u>. The dependence turns out to be only significant in 
the small vicinity of O = u>' /u> — {0, 1}, elsewhere the 
contribution being small: O (1/kl). We note that for the 
mentioned ratios fl one of the source frequencies u' or 
\oj — oj'\ must be close to zero, and the phenomenon is 
of small practical importance from experimental point 
of view. For these near-zero frequencies the wavelength 
becomes comparable to or greater than the diameter of 
the solid. Thus the substitution of the exact Green's 
function G by the Green's function in the unbounded 
medium G°° is no longer valid, and conditions for the 
time scales made in Section III are not met. 

With these limitations in mind we set X/l — for 
practical calculations, and obtain analytical expression 
for IIA3JI . As the result of this procedure a singular- 
ity of the function N at lj'/ui = {0,1} is acquired (not 
shown in Figure . Discontinuities of the function at 
the wavenumber resonances, and discontinuities of its 
slope at lo/lu' = 2/(1 + q/q) and 1 — 2/ (1 + q/c 4 ), are 
found as well. We expect in practice to observe sharp 
transitions over finite ranges in lu/lu' of order 1/kl, at a 
characteristic frequency ratio oj/uj' possibly shifted by an 
amount O {1/kl). 



